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1.  Introduction 


Energetic  molecular  crystals  are  found  in  military  and  industrial  applications  involving  explosives 
and  propellants.  For  example,  plastically  bonded  explosives  consist  of  energetic  crystals 
embedded  in  a  waxy  binder  phase  (4).  Defects  in  energetic  materials  are  thought  to  affect  their 
initiation  sensitivity.  Local  stresses  may  concentrate  in  the  vicinity  of  cracks,  pores,  or  lattice 
defects,  including  dislocations  and  point  defects,  which  may,  in  turn,  affect  initiation  of  reactions 
associated  with  burning  or  detonation  (5).  In  anisotropic  single  crystals,  availability  of  slip 
systems  associated  with  glissile  dislocations  may  lower  peak  stresses  and  decrease  sensitivity  to 
shock  initiation  (6). 

The  focus  of  this  work  is  the  mechanical  behavior  of  the  energetic  material  cyclotrimethylene 
trinitramine  (CaH^N^Oc-,),  referred  to  as  RDX  (Research  Development  Explosive).  Single 
crystals  of  RDX  belong  to  an  orthorhombic  space  group  with  eight  molecules  per  unit  cell. 
Dislocations  in  RDX  have  been  characterized  using  etch  pit  (7)  and  x-ray  topographic  (5,  9) 
techniques.  Likely  slip  systems  for  dislocation  glide  in  RDX  have  been  suggested  by  analysis  of 
anisotropic  hardness  profiles  (10)  and  indentation  experiments  (3,  11).  The  latter  experiments 
(3,  11)  also  provide  an  estimate  of  the  critical  resolved  shear  stress  associated  with  slip  initiation, 
thought  to  be  on  the  order  of  the  theoretical  strength  (i.e.,  ~  G/10  —  G/ 20,  with  G  a 
representative  elastic  shear  modulus),  which  corresponds  physically  to  homogeneous  dislocation 
nucleation.  Inelastic  behavior  of  RDX  single  crystals  has  also  been  probed  using  shock 
experiments  (6,  12)  and  molecular  dynamics  simulations  (12-14).  RDX  undergoes  a 
polymorphic  phase  transition  (a  — >  y)  at  pressures  on  the  order  of  4  GPa  (15,  16);  this  phase 
transformation  has  also  been  studied  using  atomic  models  (17,  18). 

Continuum  crystal  plasticity  theory  permits  predictive  mesoscale  modeling  of  materials  behavior 
at  length  scales  larger  than  feasible  using  atomic  or  molecular  models,  but  with  greater  resolution 
than  that  afforded  by  macroscopic  elastic-plastic  models  that  omit  description  of  anisotropy  and 
slip  system  activity.  Grain  interactions  can  be  studied  in  direct  numerical  simulations  via  finite 
element  models,  wherein  each  crystal  of  a  polycrystal  is  resolved  geometrically.  Crystal 
plasticity  models  have  been  implemented  elsewhere  to  study  shock  loading  of  the  energetic 
materials  cyclotetramethylene  tetranitramine  (HMX)  (19)  and  pentaerythritol  tetranitrate  (PETN) 
(20,  21).  Aims  of  the  present  work  are  to  develop  and  implementent  a  crystal  plasticity  model  for 
RDX. 
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The  single  crystal  elastic-plastic  model  developed  in  the  present  work  extends  the  general 
framework  of  Becker  (22)  formulated  for  cubic  crystals  loaded  to  possibly  high  pressures.  Here, 
the  model  is  applied  to  crystals  with  orthorhombic  symmetry  characteristic  of  RDX.  Anisotropic 
elastic  constants  and  pressure-dependent  compressibility  are  obtained  from  the  experimental 
literature  ( 1 ,  2,  23).  Six  slip  systems  (from  four  different  families  of  systems)  are  implemented 
following  analysis  of  indentation  loading  profiles  and  surface  impressions  ( 3 ,  10).  The  model  is 
applied  to  study  indentation,  with  a  spherical  indenter,  of  (001),  (021),  and  (210)  faces  of  single 
crystals  of  RDX.  Results  are  compared  with  numerical  predictions  from  nonlinear  isotropic 
elasticity  and  analytical  isotropic  elastic  (i.e.,  Hertzian)  solutions.  Critical  shear  strengths 
associated  with  slip  initiation  are  assigned  upon  comparison  of  results  with  load  excursion  data 
from  experiments  (3).  Bulk  slip  system  activity  is  quantified,  and  residual  surface  slip  activity  is 
compared  with  experimental  observations. 

This  report  is  organized  as  follows.  Constitutive  theory  and  material  properties  are  described  in 
section  2.  Indentation  simulations  are  reported  in  section  3.  Conclusions  follow  in  section  4. 
Notation  of  continuum  mechanics  is  used.  Boldface  type  is  used  for  vectors  and  tensors,  which 
are  referred  to  a  fixed  Cartesian  coordinate  frame  in  the  reference  and  current  configurations. 


2.  Theory 


A  general  theory  for  mechanical  behavior  of  elastic-plastic  single  crystals  is  developed  in  section 
2.1.  Features  and  requisite  properties  particular  to  RDX  are  listed  in  section  2.2. 

2.1  Single  Crystal  Model 

Let  x  —  X  (X,  t)  denote  the  motion  of  material  points  of  the  body.  The  deformation  gradient  is 

V%  =F  =  FeFp,  (1) 

where  V(-)  denotes  the  material  gradient  (i.e.,  FaA  =  V^Xa  —  c)xa/r)XA  in  component  notation), 
FE  denotes  thermoelastic  deformation  of  the  crystal  lattice,  and  Fp  represents  plastic  deformation 
due  to  slip.  Let  a  superposed  dot  denote  the  material  time  derivative  (i.e.,  d  /dt  at  fixed  X),  and 
let  superscript  —  1  denote  inversion.  The  spatial  velocity  gradient  is 

VX(VX)-1  =FF  1  =FEFE~l  +FeLpFe~\  (2) 
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where  the  plastic  velocity  gradient  associated  with  slip  rates  yk,  reference  slip  directions  Sq,  and 
reference  slip  plane  normals  m q  on  slip  systems  k  is 

Lp  =  FPFp-]  =Y,yk4®”io-  (3) 

k 

The  tensor  product  is  ®.  Slip  system  geometry  is  pushed  forward  to  the  current  configuration  via 

sk  =  FEsl  mk=FE~Tml  (4) 

where  superscript  T  denotes  the  transpose.  From  the  orthonormality  of  slip  directions  and  slip 
plane  normals,  plastic  deformation  is  isochoric: 

detFp  =  detFptrLp  =  0,  J  =  det F  =  det FE  >  0.  (5) 

The  trace  of  a  second-order  tensor  is  denoted  by  tr(-).  By  the  polar  decomposition  theorem,  let 
the  thermoelastic  deformation  be  split  as  follows: 

Fe=ReUe1  Ue=Uet,  Rei=Ret,  det/?£  =  1.  (6) 


A  logarithmic  thermoelastic  strain  measure  is  defined  as 

e  =  In  UE  =  ^ln  (FetFe),  tre  =  ln[det(t/£)]  =  Ini.  (7) 


The  logarithmic  thermoelastic  strain  is  split  into  deviatoric  and  volumetric  parts  as 

e=e/  +  ^el,  £  =  tre  =  lni, 


(8) 


where  1  is  the  second-order  unit  tensor.  Note  that  £  =  0  when  7=1.  Let  a  denote  the  usual 
Cauchy  stress  tensor.  The  stress  tensor  in  the  defoimed  but  unrotated  crystal  coordinate  system  is 


S  =  Re  1gRe. 


(9) 


The  nonlinear  elastic  component  of  the  model  follows  (22).  Henceforward,  only  the  isothermal 
case  is  considered,  since  isothermal  conditions  are  an  appropriate  assumption  for  applications  of 
the  model  discussed  in  section  3.  Let  the  unrotated  stress  be  split  into  deviatoric  and  hydrostatic 
parts: 

S  =  S/-t-51,  S  =  ^trS  =  ^tro  =  —  p,  (10) 
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with  p  the  Cauchy  pressure.  The  following  operator  extracts  the  deviatoric  part  of  a  second-order 
tensor: 

^ABCD  =  §AC§BD  ~  ^AB§CD-  (H) 

The  constitutive  equations  for  the  deviatoric  stress  and  pressure  are  (22) 

S,  =  I':C:e,+  ^(I':C:l)e,  (12) 

P  =  :C:l.  (13) 

Here,  Bq  and  B'  are  the  reference  bulk  modulus  and  the  pressure  derivative  of  the  bulk  modulus, 
and  C  is  the  tensor  of  second-order  elastic  constants,  referred  to  the  (unrotated)  crystal  frame. 

The  colon  denotes  contraction  over  two  pairs  of  indices,  e.g.,  (C  :  e')ab  =  ^abcd^cd-  Pressure 
dependence  of  shear  elastic  coefficients  implemented  elsewhere  for  cubic  crystals  (22)  is  omitted 
in  equation  13  because  of  limited  experimental  data  for  the  material  of  present  interest. 

The  second  term  on  the  right  side  of  equation  12  accounts  for  possible  effects  of  volume  changes 
on  deviatoric  stress  in  anisotropic  crystals.  The  second  term  on  the  right  side  of  equation  13 
accounts  for  possible  effects  of  deviatoric  strains  on  pressure  in  anisotropic  crystals.  Viscoelastic 
effects  are  omitted. 


The  flow  rule  for  slip  is  (22) 


xk  =  G  :  (sk®mk) 


Tk 

T0 


Iya1 

Yo 


sgn(y* 


(14) 


The  resolved  shear  stress  on  slip  system  k  is  Tk.  Material  parameters  are  initial  and  constant  slip 
strength  Tq  for  each  slip  system,  reference  strain  rate  jo,  and  rate  sensitivity  m.  Constant  £,  -C  Yo 
provides  a  finite  strength  at  zero  strain  rate  and  is  included  for  numerical  reasons.  Strain 
hardening  is  omitted  in  equation  14  because  sufficient  data  do  not  exist  to  construct  more 
complex  slip  system  constitutive  equations  for  the  material  of  present  interest. 


2.2  RDX 


Physical  properties  of  RDX  single  crystals  are  listed  in  table  1.  This  description  is  restricted  to 
the  a  phase,  which  is  the  stable  polymorph  for  pressures  under  ^3.8  GPa  and  temperatures  under 
m  480  K.  High  pressure  phases  with  different  structures  (y,  |3 )  and  melting  at  higher  temperatures 
are  not  considered  relevant  to  this  application  and  are  discussed  elsewhere  (75). 
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Table  1.  Structural  and  physical  properties  of  RDX  single  crystals  (ambient) . 


Property 

Value 

Reference 

Space  group 

Pbca 

(10) 

Crystal  structure 

orthorhombic 

Lattice  parameters  [nm] 

a=  1.3182 

b=  1.1574 

c=  1.0709 

(10) 

Mass  density  [g/cnr] 

1.806 

(2) 

Elastic  properties  are  listed  in  table  2,  obtained  from  experiments.  Isentropic  second-order 
elastic  constants  ( 1 ,  2)  and  the  bulk  modulus  have  been  converted  to  isothermal  values  at  295  K 
via  the  usual  thermodynamic  formulae  (24)  incorporating  anisotropic  thermal  expansion  (25)  and 
specific  heat  (17).  Voigt’s  notation  is  used,  i.e.,  C abcd  ^  Cap ,  where  Greek  indices  1 , 2, . . . ,  6. 
Voigt  (Gy)  and  Reuss  (Gr)  bounds  (24,  26)  on  the  effective  shear  modulus  are  also  listed. 
Differences  between  Voigt  and  Reuss  bounds  for  the  bulk  modulus  Bq,  on  the  order  1  —  3%,  are 
considered  insignificant. 

Table  2.  Isothermal  second-order  elastic  constants  of  RDX  single  crystals  (ambi¬ 
ent). 


Property 

Value  (1) 

Value  (2) 

Cn  [GPa] 

24.56 

36.48 

C22  [GPa] 

18.85 

24.49 

C33  [GPa] 

17.33 

20.78 

Cn  [GPa] 

7.61 

0.90 

C13  [GPa] 

5.30 

1.26 

C23  [GPa] 

5.24 

8.16 

C44  [GPa] 

5.15 

11.99 

C55  [GPa] 

4.06 

2.72 

C66  [GPa] 

6.90 

7.68 

Bq  [GPa] 

10.5 

11.2 

Gv  [GPa] 

6.06 

9.26 

Gr  [GPa] 

5.72 

6.40 

As  is  evident  from  values  listed  in  table  2,  second-order  elastic  constants  reported  for  different 
sets  of  experiments  on  single  crystals  of  RDX  can  vary  substantially.  Values  obtained  using 
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resonant  ultrasonic  methods  listed  in  table  2  (i)  are  in  reasonably  close  agreement  with  those 
reported  by  other  researchers  using  the  same  method  (27).  Values  obtained  using  Brillouin 
scattering  (2)  listed  in  table  2  are  notably  different,  with  particularly  larger  shear  stiffness.  Values 
obtained  using  a  third  technique,  impulsive  stimulated  thermal  scattering  (28),  are  similar  to  those 
obtained  using  resonant  ultrasound  spectroscopy.  Values  obtained  using  empirical  atomic  models 
(18,  29)  also  display  differences  among  one  another  and  differences  with  those  obtained  in 
experiments. 

In  this  work,  the  two  sets  of  elastic  constants  of  Haussuhl  (1)  and  Haycraft  et  al.  (2)  are  used 
because  these  are  the  softest  and  stiffest  reported  experimental  measurements,  respectively; 
results  obtained  using  these  two  sets  are  thought  to  provide  an  indication  of  the  most  and  least 
compliant  expected  response.  Results  obtained  using  the  constants  of  Haussuhl  (1)  are  very 
similar  to  those  that  would  be  obtained  using  values  from  Schwarz  et  al.  (27)  and  Sun  et  al.  (28). 

Pressure  and  temperature  dependencies  of  second-order  elastic  coefficients  have  been  calculated 
using  molecular  dynamics  (29);  however,  these  predicted  values  have  not  been  validated  using 
experiments,  and  some  discrepancies  exist  among  calculated  second-order  elastic  constants  at 
room  temperature  and  experimental  values  (2).  For  this  reason,  in  this  study  the  nonlinear  elastic 
model  only  incorporates  pressure  dependence  of  the  compressibility  (i.e.,  the  bulk  modulus) 
obtained  experimentally  and  not  that  of  all  elastic  coefficients.  The  value  used  here  is  B'  —  6.95 
(23,  30). 

Potential  slip  systems  in  RDX-as  identified  from  hardness  versus  orientation  profiles  (10), 
indentation  force  versus  depth  data  (3,  11),  and  residual  surface  impressions  from  indentation  (3)- 
are  listed  in  table  3.  Slip  system  geometry  is  referred  to  a  Cartesian  system  with  axes  (Xi,X2,Xi) 
parallel  to  lattice  vectors  (a,b,c).  These  slip  systems  are  illustrated  in  figure  1. 


Table  3.  Slip  systems  in  RDX  single  crystals  (indentation). 


System  k 

Miller  indices 

m0 

•so 

Maximum  strength  [GPa] 

Reference 

1 

(021)  [100] 

(0,0.880,0.475) 

[1,0,0] 

0.885 

(3,  10) 

2 

(021)  [100] 

(0,-0.880,0.475) 

[1,0,0] 

3 

(011)  [100] 

(0,0.679,0.734) 

[1,0,0] 

0.645 

(3,  10) 

4 

(Oil)  [100] 

(0,-0.679,0.734) 

[1,0,0] 

5 

(010)  [100] 

(0,1,0) 

[1,0,0] 

0.885 

(3) 

6 

(010)  [001] 

(0,1,0) 

[0,0,1] 

0.885 

(3,  10) 

Listed  initial  slip  system  strengths  are  upper  bounds  obtained  from  analysis  of  load  excursion  data 
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{021 }  [100]  slip 
2  planes  x  1  direction 


(010)[100]  slip 
1  plane  x  1  direction 


X,[100] 


(010)  [001]  slip 
1  plane  x  1  direction 


Figure  1.  Slip  systems  in  RDX  single  crystals  (unit  cell  parameters  not  to  scale). 


using  the  analytical  Hertzian  solution  for  frictionless  spherical  indentation  into  a  semi-infinite, 
linear  elastic,  isotropic  material.  It  is  noted  that  other  systems  may  become  active  (and  those 
listed  may  become  inactive)  for  loading  regimes  involving  very  different  pressures,  temperatures, 
and/or  strain  rates.  For  example,  molecular  dynamics  simulations  (12)  suggest  that  partial 
dislocation  loops  may  glide  on  (001)[010]  during  shock  loading  to  pressures  in  excess  of  ~  1  GPa. 

In  numerical  simulations  of  section  3,  shear  strengths  for  various  families  of  systems  are  varied 
parametrically  between  physically  reasonable  bounds  on  the  order  of  the  theoretical  strength  as 
suggested  in  the  literature  (3,  11),  where  G  —  Gr  —  6.4  GPa  (2).  These  and  other  parameters  are 
listed  in  table  4.  Reference  strain  rate  jo  is  on  the  order  of  the  plastic  slip  rate  experienced  during 
indentation  simulations  described  later,  and  the  very  small  value  of  rate  sensitivity  m  provides  for 
a  nearly  rate  insensitive  response,  the  simplest  numerically  pragmatic  assumption  in  the  absence 
of  experimental  data  enabling  calibration  of  this  material  parameter. 


Table  4.  Crystal  plasticity  model  parameters  for  RDX. 


Property 

Value 

xo 

G/40  <  Tq  <  G/10 

m 

5  x  1(T3 

jo 

<N 

1 

o 

% 

o 

1 
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3.  Indentation 


The  indentation  boundary  value  problem  and  numerical  aspects  are  described  in  section  3.1. 
Results  for  indentation  of  RDX  with  slip  suppressed  (i.e.,  nonlinear  anisotropic  elasticity)  are 
reported  in  section  3.2.  Results  for  indentation  of  RDX  with  slip  enabled  are  reported  in  section 
3.3. 


3.1  Boundary  Value  Problem 

The  constitutive  model  presented  in  section  2  is  implemented  in  the  ALE  3D  multi -physics  code. 
Simulations  of  indentation  are  performed  using  an  implicit  solver  for  static  equilibrium  (i.e.,  no 
inertial  terms  in  the  balance  of  linear  momentum). 

The  problem  geometry  is  intended  to  mimic  previous  experimental  studies  (3,  11).  A  spherical 
diamond  indenter  of  radius  R,  with  a  nominal  value  of  R=  1.482  (im,  is  used  to  indent  a  flat 
surface  of  a  single  crystal  of  RDX  of  variable  lattice  orientation.  Diamond  is  represented  as  an 
isotropic  nonlinear  elastic  material  with  properties  listed  in  table  5. 


Table  5.  Isotropic  elastic  properties  of  diamond  indenter. 


Property 

Value 

Reference 

B0  [GPa] 

443 

(31) 

Gv  [GPa] 

538 

B' 

4.0 

(32) 

Following  Nowak  et  al.  (33),  the  substrate  is  represented  by  a  right  circular  cylinder.  The 
cylinder  is  assigned  a  height  and  radius  of  2 R;  further  increases  in  dimensions  of  the  cylinder  did 
not  affect  results  of  interest.  The  indenter  is  modeled  as  a  half-sphere.  Each  body  (i.e.,  indenter 
and  substrate)  is  discretized  using  eight-node  hexahedral  elements  with  selective  reduced 
integration;  the  entire  problem  consists  of  ~  55000  finite  elements.  The  mesh  of  the  substrate  is 
highly  refined  in  the  vicinity  of  contact  beneath  the  indenter,  where  stress  fields  are 
inhomogeneous,  and  the  mesh  coarsens  progressively  with  distance  from  the  initial  contact  point. 
Simulations  with  smaller  elements  demonstrated  that  further  increases  in  mesh  refinement  did  not 
affect  results  of  interest.  Anisotropic  crystal  plasticity  simulations  were  also  performed  with  the 
original  cylindrical  mesh  rotated  by  n/4  about  the  X3  axis  (i.e.,  Z-axis)  with  the  initial  crystal 
lattice  orientation  held  fixed;  nearly  identical  results  were  obtained  for  plastic  slip  contours, 


8 


demonstrating  that  strain  localization  was  not  artificially  sensitive  to  mesh  construction. 

During  the  loading  phase,  the  upper  face  of  the  half-spherical  indenter  is  assigned  a  constant 
(downward)  velocity  of  D  =  10  nm/s,  leading  to  strain  rates  similar  to  those  of  experiments  (3). 
The  lower  face  of  the  cylinder  is  rigidly  fixed,  while  the  lateral  sides  (circumference)  are  traction 
free.  The  indentation  depth  is  denoted  by  D;  actual  depth  d  of  the  tip  of  the  sphere  in  contact 
with  the  surface  is  monitored  as  an  outcome  of  the  numerical  solution  procedure.  Only  for  a  rigid 
indenter  would  d  —  D.  Indentation  force  P  is  defined  as  the  sum  of  nodal  forces  along  the  upper 
face  of  the  half  sphere  acting  in  the  direction  of  D,  i.e.,  the  sum  of  forces  work  conjugate  to  the 
prescribed  nodal  velocities. 

In  most  of  the  simulations,  contact  between  indenter  and  substrate  is  assumed  frictionless, 
following  previous  studies  that  rely  on  analytical  solutions  (3,  11).  In  some  simulations  discussed 
later,  unloading  is  also  performed,  whereby  after  a  peak  depth  is  attained,  the  upper  face  of  the 
indenter  is  assigned  an  upward  velocity  of  D  —  10  nm/s  until  contact  is  released.  If  plastic 
deformation  has  occurred,  some  residual  deformation  remains  in  the  substrate  upon  unloading. 

Comparison  of  simulation  results  (and  published  experimental  results)  with  the  available 
analytical  solution  proves  useful.  The  Hertzian  solution  for  spherical  indentation  of  a  linear 
elastic  isotropic  half  space  is  (3) 

P=^ErR1/2d3/2,  (15) 

where  the  reduced  elastic  modulus  is 

Er=  [(l-vft/E.  +  tl-v^/Ej-1,  (16) 

with  v,  ( V.v )  and  E,  (Es)  denoting  Poisson’s  ratio  and  Young’s  modulus,  respectively,  of  the 
indenter  (substrate).  In  the  present  context,  these  are  computed  from  Bq  and  G  using  the  standard 
formulae  (24).  The  solution  is  fairly  insensitive  to  properties  of  diamond  because  its  stiffness 
exceeds  that  of  RDX  by  1-2  orders  of  magnitude.  It  is  noted  that  equation  15  neglects  nonlinear 
aspects  of  the  solution  associated  with  large  changes  in  geometry  and  nonlinear  material  behavior 
(e.g.,  finite  volume  and  shape  changes  and  nonlinear  elastic  material  properties).  Effects  of 
anisotropy  and  plasticity  are,  of  course,  also  omitted  in  the  analytical  solution. 

Simulations  of  indentation  onto  (001),  (021),  and  (210)  planes  are  reported  in  the  following. 
Rotation  matrices  and  second-order  elastic  constants  associated  with  these  crystal  orientations  are 
listed  in  the  appendix. 
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3.2  Elastic  Indentation 


Isotropic  elastic  simulations  are  considered  first  in  order  to  facilitate  comparison  with  the 
analytical  solution,  equation  15.  Figures  2  and  3  show  indentation  force,  P,  versus  tip  depth,  d, 
for  indentation  simulations  incorporating  nonlinear  isotropic  elasticity  with  Voigt  (Gy)  and  Reuss 
( Gr )  bounds  on  the  shear  modulus.  Also  shown  are  predictions  using  Gr  with  a  constant  bulk 
modulus  (B'  =  0).  Results  in  figure  2  correspond  to  elastic  constants  of  Haussuhl  ( 1 );  results  in 
figure  3  correspond  to  elastic  constants  of  Haycraft  et  al.  (2).  The  indenter  radius  is  R  —  1.482 
pm  (3).  Analytical  solutions  using  the  same  bounds  on  the  shear  modulus,  and  incorporating 
linear  elastic  behavior  as  in  equation  15,  are  also  plotted.  As  expected,  force  is  larger  when  the 
larger  shear  modulus  (Gy)  is  used.  For  tip  depth  d  >50  nm,  the  pressure  dependence  of  the  bulk 
modulus  causes  an  increase  in  force  relative  to  the  case  with  B'  —  0. 


Figure  2.  Indentation  force  vs.  tip  depth:  isotropic  elastic  simulation  results  and 
analytical  solutions,  R  =  1 .482  pm,  elastic  constants  of  Haussuhl  (1). 
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Figure  3.  Indentation  force  vs.  tip  depth:  isotropic  elastic  simulation  results  and 
analytical  solutions,  R  =  1.482  pm,  elastic  constants  of  Haycraft  et 
al.  (2). 


Differences  between  analytical  and  numerical  solutions  become  noticeable  for  d  >  50  nm  when 
Gr  is  used,  and  for  d>  40  nm  when  Gy  is  used.  In  both  cases,  numerical  predictions  are  stiffer 
than  corresponding  analytical  solutions.  Differences  between  numerical  and  analytical  solutions 
are  attributed  to  the  following  features  resolved  by  the  calculation  that  are  not  included  in  the 
analytical  solution:  (1)  geometric  nonlinearities  associated  with  problem  geometry,  (2)  finite 
logarithmic  deviatoric  strain  measures  and  true  volume  changes,  (3)  nonlinear  compressibility  for 
cases  with  B'  >  0,  and  (4)  finite  size  of  the  substrate.  Recall  from  discussion  in  section  3.1  that 
mesh  construction  did  not  influence  results  such  as  indentation  force  P;  differences  between 
numerical  and  analytical  solutions  are  not  attributed  to  mesh  resolution.  All  results  discussed 
hereafter  involve  nonlinear  compressibility. 


Closer  agreement  between  the  numerical  model  (with  constant  bulk  modulus)  and  linear  elastic 
analytical  solutions  can  be  obtained  by  the  transformation  d  — >  d  +  8  in  the  analytical  solution  of 
equation  15,  where  8  is  a  correction  that  accounts  for  the  presence  of  the  rigid  boundary  at  a 
depth  Z  —  2R  from  the  surface  of  the  substrate: 


8  = 


P 

2nZEs 


(3  +  vs 


(17) 


As  shown  in  figure  4,  equation  17  yields  reasonably  close  agreement  between  model  and 
(corrected)  analytical  solutions,  suggesting  that  differences  between  numerical  predictions  and 
uncorrected  analytical  solutions  may  be  due,  at  least  in  part,  to  the  finite  size  of  the  substrate 
modeled  in  the  simulations.  It  should  be  noted  that  this  correction,  obtained  from  the  analytical 
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Boussinesq  solution  for  displacement  due  to  a  point  force  applied  to  a  linear  elastic  isotropic 
half-space  (34),  is  approximate. 


Figure  4.  Indentation  force  vs.  tip  depth:  isotropic  elastic  simulation  results  and 
analytical  solutions  with  and  without  corrections,  R  =  1.482  pm;  soft 
elastic  constants  of  Haussuhl  (i)  and  stiff  constants  of  Haycraft  et  al.  (2), 
Reuss  shear  modulus. 


Results  from  simulations  incorporating  isotropic  and  anisotropic  nonlinear  elasticity  are 
compared  in  figures  5  and  6,  which  show  force  P  versus  applied  indentation  depth  D.  The 
indenter  radius  is  R  =  1.482  pm  (3).  The  difference  between  tip  depth  d  and  applied 
displacement  D  is  ~  1-2%  over  the  duration  of  each  simulation.  Load  excursion  data  thought  to 
correspond  to  the  onset  of  slip  in  indentation  experiments  (3)  are  also  shown. 
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D  [nm] 


150 


Figure  5.  Indentation  force  vs.  applied  displacement:  isotropic  and  anisotropic 
elastic  simulation  results,  and  experimental  initial  load  excursion  data 
(3),  R  =  1.482  pm,  elastic  constants  of  Haussuhl  (i). 


D  [nm] 

Figure  6.  Indentation  force  vs.  applied  displacement:  isotropic  and  anisotropic 
elastic  simulation  results,  and  experimental  initial  load  excursion  data 
(, 3 ),  R  =  1 .482  pm,  elastic  constants  of  Haycraft  et  al.  (2). 
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First,  consider  the  numerical  predictions  in  figure  5,  which  are  obtained  using  the  relatively 
compliant  elastic  constants  reported  by  Haussuhl  (/).  Orientation  (210)  is  relatively  stiff,  with 
predicted  force  P  exceeding  that  of  the  isotropic  case  using  Voigt  shear  modulus  Gy.  Orientation 
(001)  is  the  most  compliant,  with  predicted  force  P  less  than  that  for  the  isotropic  case  using 
Reuss  shear  modulus  Gr.  Simulations  predict  p(°01)  <  p(021)  <  p(210).  As  shown  in  the 
appendix,  the  elastic  stiffness  component  associated  with  indentation  follows  the  sequence 
eg01)  <  eg21)  <  eg10).  Experiments  (5)  suggest  that  orientation  (001)  is  most  compliant,  in 
agreement  with  the  present  results.  Experiments  also  suggest  that  (210)  has  a  higher  effective 
modulus  Er  than  (021),  while  (021)  has  a  higher  measured  hardness  than  (210).  Now  compare 
numerical  predictions  with  experimental  data  in  figure  5.  Close  agreement  between  model 
predictions  and  load  excursion  data,  within  the  accuracy  of  the  data,  e.g.,  «  1 0-20%  uncertainty  in 
measured  P  (3),  is  demonstrated  for  all  three  orientations  of  the  RDX  substrate. 

Next  consider  the  numerical  predictions  in  figure  6,  which  are  obtained  using  the  relatively  stiff 
elastic  constants  reported  by  Haycraft  et  al.  (2).  Orientations  (021)  and  (210)  are  relatively  stiff, 
with  predicted  force  P  comparable  to  the  isotropic  case  using  Voigt  shear  modulus  Gy. 
Orientation  (001)  is  relatively  compliant,  with  predicted  force  P  slightly  exceeding  the  isotropic 
case  using  Reuss  shear  modulus  Gr.  In  other  words,  simulations  predict  p(°oi)  <-  p(02i)  ^  p(2i0) 
Next  compare  numerical  predictions  with  experimental  data  in  figure  6.  For  the  (210) 
orientation,  the  load  curve  from  the  elastic  model  coincides  with  the  load  excursion  data,  within 
the  accuracy  of  the  data.  Elastic  model  load  curves  are  significantly  stiffer  than  load  excursion 
data  for  (001)  and  (021)  orientations.  Molecular  simulations  suggest  that  all  shear  moduli  of 
RDX  increase  with  increasing  pressure  (29);  incorporation  of  such  effects  into  the  present 
indentation  simulations  would  seem  to  lead  to  an  increase  in  tangent  modulus  with  increasing 
depth  of  indentation,  further  increasing  differences  between  modeling  using  constants  of  Haycraft 
et  al.  (2)  and  experiment. 

Predictions  obtained  using  each  set  of  elastic  constants  (1,  2)  are  compared  directly  for  substrate 
orientations  (001),  (021),  and  (210)  in  figures  7,  8,  and  9.  Significantly  closer  agreement  between 
model  and  experiment  (2)  is  obtained  from  the  elastic  constants  of  Haussuhl  ( 1 )  for  indentation 
onto  (001)  and  (021)  planes.  Comparable  accuracy  is  obtained  from  either  set  of  elastic  constants 
for  indentation  on  (210)  planes.  Thus,  the  present  simulations  confirm  that  elastic  constants 
obtained  using  resonant  ultrasonic  methods  ( 1 ,  27)  provide  a  more  realistic  representation  of 
elastic  stiffness  during  nano-indentation  than  elastic  constants  obtained  from  Brillouin  scattering 
(2),  with  the  latter  appearing  too  stiff. 
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2.5 


Figure  7.  Indentation  force  vs.  applied  displacement  for  elastic  indentation  onto 
plane  (001):  predictions  obtained  using  stiff  (2)  and  soft  (i)  elastic  con¬ 
stants  compared  to  experiment  ( 3 ). 


D  [nm] 

Figure  8.  Indentation  force  vs.  applied  displacement  for  elastic  indentation  onto 
plane  (021):  predictions  obtained  using  stiff  (2)  and  soft  (i)  elastic  con¬ 
stants  compared  to  experiment  ( 3 ). 
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D  [nm] 

Figure  9.  Indentation  force  vs.  applied  displacement  for  elastic  indentation  onto 
plane  (210):  predictions  obtained  using  stiff  (2)  and  soft  ( 1 )  elastic  con¬ 
stants  compared  to  experiment  ( 3 ). 
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The  influence  of  indenter  radius  R  on  nonlinear  anisotropic  simulation  results  is  shown  in  figures 
10,  11,  and  12;  in  each  case,  elastic  constants  of  Haycraft  et  al.  (2)  are  employed.  For  all  three 
orientations  of  the  RDX  substrate,  indentation  force  is  nearly  proportional  to  R1/2,  as  in  the 
isotropic  linear  elastic  analytical  solution  equation  15.  For  exact  R1'2  proportionality,  the  three 
curves  within  each  of  figures  10,  11,  and  12  would  coincide.  This  supports  the  assertion  that 
nonlinearities  contribute  to  differences  between  numerical  and  analytical  solutions. 


Figure  10.  Normalized  force  vs.  tip  depth  for  elastic  indentation  onto  plane  (001), 
elastic  constants  of  Haycraft  et  al.  (2). 


Figure  11.  Normalized  force  vs.  tip  depth  for  elastic  indentation  onto  plane  (021), 
elastic  constants  of  Haycraft  et  al.  (2). 
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Figure  12.  Normalized  force  vs.  tip  depth  for  elastic  indentation  onto  plane  (210), 
elastic  constants  of  Haycraft  et  al.  (2). 


All  simulation  results  reported  in  figures  in  this  work  are  obtained  assuming  frictionless  contact. 
In  order  to  estimate  an  upper  bound  on  indentation  force  for  non-negligible  sliding  friction 
between  indenter  and  substrate,  simulations  incorporating  nonlinear  isotropic  elasticity  with  Gy 
or  Gr  of  Haycraft  et  al.  (2)  were  also  conducted  assuming  fully  tied  tangential  contact,  i.e.,  no 
sliding,  full  sticking,  or  an  infinite  friction  coefficient.  In  such  simulations,  normal  separation 
was  permitted.  It  was  found  that  differences  between  frictionless  and  tied  contact  were 
insignificant.  For  example,  at  D  =  150  nm,  for  simulations  incorporating  Gy,  P  =  2.56  mN  for 
frictionless  contact  and  P  =  2.57  mN  for  tied  contact.  Experimental  measurements  of  dynamic 
friction  for  RDX  single  crystals  sliding  on  a  glass  substrate  (35)  suggest  a  friction  coefficient  on 
the  order  of  unity  for  loads  under  1  g  (~  10  mN)  wherein  contact  is  characterized  as  elastic,  and  a 
friction  coefficient  of  0.35  for  loads  in  excess  of  10  g,  wherein  contact  is  characterized  as  plastic. 
Experimental  sliding  velocities  were  on  the  order  of  0.2  mm/s.  The  present  indentation 
simulations  consider  a  different  geometry,  smaller  system  sizes  (loads  under  1  g),  and  slower 
sliding  velocities  (on  the  order  of  10  nm/s),  so  the  reported  experimental  values  for  friction 
coefficients  (35)  may  not  strictly  apply  here. 

3.3  Elastic-Plastic  Indentation 

Results  are  now  reported  for  indentation  of  RDX  single  crystals  described  by  the  complete 
nonlinear  crystal  plasticity  constitutive  model  of  section  2.  Because  strengths  of  families  of  slip 
systems  are  not  known  precisely  a  priori  from  experiments,  shear  strengths  Tq  are  varied 
parametrically  over  the  range  listed  in  table  4.  This  range,  which  corresponds  to  slip  resistance 
on  the  order  of  the  theoretical  strength,  is  physically  descriptive  of  homogeneous  nucleation  of 
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glissile  dislocation  lines  and  loops. 


Results  of  indentation  force  versus  applied  indention  depth  are  shown  in  figures  13,  14,  and  15. 

In  each  of  these  figures,  results  correspond  to  elastic  constants  of  Haycraft  et  al.  (2).  For 
situations  in  which  one  family  (consisting  of  one  or  two  systems,  as  in  figure  1)  is  assigned 
strength  x/j  =  G/ 20,  the  remaining  systems  are  assigned  strength  G/10.  Recall  that  Miller 
indices  corresponding  to  slip  system  numbers  1-6  are  listed  in  table  3.  For  indentation  on  (001), 
only  slip  system  6  exhibits  significant  activity;  lowering  the  strength  of  any  of  the  other  systems 
1-5  from  G/ 10  to  Gj 20  does  not  noticeably  affect  indentation  force  P.  For  indentation  on  (021), 
lowering  the  strength  of  system  6  to  G/20  provides  the  greatest  reduction  in  indentation  force. 
However,  lowering  the  strength  of  the  other  families  of  systems  does  produce  a  slight  decrease  in 
force.  For  indentation  on  (210),  the  opposite  trend  is  observed.  Lowering  the  strength  of  system 
6  to  G/20  produces  a  small  reduction  in  force,  whereas  lowering  the  strength  of  the  other  families 
of  systems  produces  a  more  substantial  reduction. 


D  [nm] 


Figure  13.  Force  vs.  depth  for  elastic -plastic  indentation  onto  plane  (001);  slip  sys¬ 
tems  have  strength  G/10  unless  noted  otherwise,  elastic  constants  of 
Haycraft  et  al.  (2). 
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Figure  14.  Force  vs.  depth  for  elastic-plastic  indentation  onto  plane  (021);  slip  sys¬ 
tems  have  strength  G  j  1 0  unless  noted  otherwise,  elastic  constants  of 
Haycraft  et  al.  (2). 


D  [nm] 


Figure  15.  Force  vs.  depth  for  elastic-plastic  indentation  onto  plane  (210);  slip  sys¬ 
tems  have  strength  G / 1 0  unless  noted  otherwise,  elastic  constants  of 
Haycraft  et  al.  (2). 
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Model  predictions  for  uniform  slip  strength  Tq  —  G/ 20  are  compared  with  experimental  data  of 
Ramos  et  al.  (3)  in  figures  16,  17,  and  18,  with  predictions  obtained  using  elastic  constants  of 
references  (7)  and  (2).  Predicted  forces  exceed  experimental  values  at  larger  indendation  depths 
in  each  orientation  and  for  both  sets  of  elastic  constants.  Results  corresponding  to  more 
compliant  elastic  constants  (/)  provide  closer  agreement  to  experimental  values  than  results 
corresponding  to  stiffer  elastic  constants  (2). 


D  [nm] 


Figure  16.  Force  vs.  depth  for  elastic-plastic  indentation  onto  plane  (001).  Model 
results  assume  all  slip  systems  have  same  strength  G/ 20;  predictions 
obtained  using  stiff  (2)  and  soft  (i)  elastic  constants;  experimental  data 
( 3 )  shown  for  comparison. 


21 


Figure  17.  Force  vs.  depth  for  elastic-plastic  indentation  onto  plane  (021).  Model 
results  assume  all  slip  systems  have  same  strength  G/ 20;  predictions 
obtained  using  stiff  (2)  and  soft  (i)  elastic  constants;  experimental  data 
(. 3 )  shown  for  comparison. 


Figure  18.  Force  vs.  depth  for  elastic-plastic  indentation  onto  plane  (210).  Model 
results  assume  all  slip  systems  have  same  strength  G/ 20;  predictions 
obtained  using  stiff  (2)  and  soft  (i)  elastic  constants;  experimental  data 
(. 3 )  shown  for  comparison. 
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Model  predictions  for  uniform  slip  strengths  Tq  of  G/10,  G/ 20,  and  G/ 40  are  compared  with 
experimental  data  (3)  in  figures  19,  20,  and  21,  with  predictions  obtained  using  elastic  constants 
of  Haycraft  et  al.  (2).  Further  reduction  of  slip  strengths  below  G/40  enables  closer  agreement 
with  experimental  data  at  larger  indentation  depths,  but  also  leads  to  significant  under-prediction 
of  the  force  at  smaller  depths  near  initial  yield. 


D  [nm] 


Figure  19.  Force  vs.  depth  for  elastic-plastic  indentation  onto  plane  (001),  elastic 
constants  of  Flaycraft  et  al.  (2).  Model  results  assume  all  slip  systems 
have  same  strength;  experimental  data  (3)  shown  for  comparison. 
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Figure  20.  Force  vs.  depth  for  elastic-plastic  indentation  onto  plane  (021),  elastic 
constants  of  Haycraft  et  al.  (2).  Model  results  assume  all  slip  systems 
have  same  strength;  experimental  data  ( 3 )  shown  for  comparison. 


D  [nm] 

Figure  21.  Force  vs.  depth  for  elastic-plastic  indentation  onto  plane  (210),  elastic 
constants  of  Haycraft  et  al.  (2).  Model  results  assume  all  slip  systems 
have  same  strength;  experimental  data  ( 3 )  shown  for  comparison. 
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As  noted  in  section  3.2,  uncertainty  in  the  true  tip  radius  R  of  the  indenter  could  lead  to 
discrepancies  between  simulations  and  experiments.  Also  noted  in  section  3.2,  surface  fractures 
and/or  subsurface  fractures  could  contribute  to  a  loss  of  stiffness  that  would  be  reflected  only  in 
the  experimental  data.  RDX  is  prone  to  cleavage  fracture  on  planes  (001),  (010),  (001),  (241), 
and  (241)  ( 3 ,  7,  9).  Experimental  data  demonstrate  nearly  horizontal  steps  in  force  versus 
displacement  corresponding  to  discrete  slip  (and  possibly  fracture)  events  that  are  not  readily 
resolved  by  a  constant  strength  continuum  crystal  plasticity  model  such  as  the  one  formulated 
here. 

Details  of  slip  system  interactions  and  pressure  dependence  of  shear  strength  are  omitted  in  the 
present  model.  Indentation  experiments  have  suggested  the  importance  of  cross  slip  (3),  and 
atomic  modeling  has  noted  that  different  slip  mechanisms  may  become  important  at  high 
pressures  (12).  Incorporation  of  these  effects  into  a  more  complex  slip  model  might  provide 
closer  agreement  with  experiment,  e.g.,  if  glide  resistance  were  to  decrease  with  pressure,  since 
local  pressures  under  the  indenter  achieve  several  GPa,  as  will  be  shown  later.  Atomic  modeling 
(12,  36)  may  provide  insight  into  dependence  of  slip  resistance  on  pressure  (and  temperature,  etc.) 
not  available  from  experimental  methods. 

In  subsequent  figures  presented  in  this  work,  elastic  constants  obtained  from  resonant  ultrasound 
spectroscopy  (1)  are  used,  and  all  systems  are  assigned  a  strength  of  x/,  =  G/ 20  =  0.32  GPa.  The 
intent  is  not  necessarily  reproduction  of  experimental  force-displacement  data,  but  rather  to 
deduce  activity  of  which  slip  system(s)  is  most  prominent  for  indentation  on  particular  planes. 

The  present  crystal  plasticity  simulations  enable  direct  quantification  of  surface  and  subsurface 
slip  on  each  system.  In  contrast,  hardness  or  indentation  experiments  (3,  10,  11)  require 
significant  interpretation  of  data  to  deduce  slip  activity  and  do  not  provide  a  quantitative  measure 
of  the  relative  contributions  of  each  slip  system  to  the  overall  plastic  strain  field.  In  the 
aforementioned  experiments,  visual  observations  of  slip  traces  are  restricted  to  residual  surface 
profiles,  whereas  simulations  enable  visualization  of  subsurface  slip  activity. 

Simulations  demonstrated  sufficient  mesh  refinement  for  elastic-plastic  simulations.  In 
particular,  doubling  of  the  number  of  elements  in  the  mesh  resulted  in  a  less  than  1%  change  in  P 
for  indentation  depths  up  to  300  nm  for  all  three  orientations  of  the  substrate  with  Xq  =  G/20 
assigned  for  all  k. 

Model  predictions  of  cumulative  slip  for  indentation  to  a  depth  of  D  —  200  nm  onto  planes  (001), 
(021),  and  (210)  are  shown  in  figures  22,  23,  and  24.  Corresponding  pressure  contours  are  shown 
in  figures  25,  26,  and  27.  In  each  case,  a  slice  along  the  centerline  of  the  cylinder  normal  to  the 
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laboratory  X\  axis  (i.e.,  X-axis)  is  shown.  For  indentation  on  (001)  and  (021)  planes  the  X\  axis 
is  normal  to  a  (100)  plane;  for  indentation  on  (210)  the  X\  axis  is  approximately  normal  to  a  (230) 
plane.  Cumulative  total  slip  y  is  defined  as 

Y  =  EYt  =  E/lV‘l<".  (18) 

k  k  J 

with  yk  the  monotonically  increasing  cumulative  slip  on  system  k.  Slip  activity  is  greater  for 
indentation  into  (021)  and  (210)  planes  than  for  indentation  on  (001).  For  this  particular  viewing 
plane,  slip  and  pressure  contours  for  indentation  on  (021)  are  noticeably  asymmetric.  Maximum 
local  pressure  (but  not  average  pressure)  is  highest  for  indentation  on  (001)  despite  the  lower 
elastic  stiffness  for  this  orientation.  The  wireframe  mesh  of  the  indenter  is  drawn  in  each  figure; 
in  order  to  enable  clear  visualization  of  heterogeneous  slip  and  pressure  distributions  in  the  RDX, 
the  mesh  of  the  substrate,  which  is  considerably  more  refined  than  that  of  the  indenter,  is  not 
shown. 
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Figure  22.  Cumulative  total  slip  contours  for  indentation  to  depth  D  =  200  nm 
onto  plane  (001);  a  slice  along  the  centerline  of  the  cylinder  normal  to 
the  laboratory  X\  axis  [i.e.,  X-axis  or  (100)]  is  shown. 
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Figure  23.  Cumulative  total  slip  contours  for  indentation  to  depth  D  =  200  nm 
onto  plane  (021);  a  slice  along  the  centerline  of  the  cylinder  normal  to 
the  laboratory  X\  axis  [i.e.,  X-axis  or  (100)]  is  shown. 
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Figure  24.  Cumulative  total  slip  contours  for  indentation  to  depth  D  =  200  nm 
onto  plane  (210);  a  slice  along  the  centerline  of  the  cylinder  normal  to 
the  laboratory  X\  axis  [i.e.,  X-axis  or  (230)]  is  shown. 


27 


pressure  (MPa) 


-2.0  -1.0  0.0  1.0  2.0 
r  juii 


Figure  25.  Hydrostatic  pressure  contours  for  indentation  to  depth  D  =  200  nm  onto 
plane  (001);  a  slice  along  the  centerline  of  the  cylinder  normal  to  the 
laboratory  X\  axis  [i.e.,  X-axis  or  (100)]  is  shown. 
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Figure  26.  Hydrostatic  pressure  contours  for  indentation  to  depth  D  =  200  nm  onto 
plane  (021);  a  slice  along  the  centerline  of  the  cylinder  normal  to  the 
laboratory  X\  axis  [i.e.,  X-axis  or  (100)]  is  shown. 
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Figure  27.  Hydrostatic  pressure  contours  for  indentation  to  depth  D  =  200  nm  onto 
plane  (210);  a  slice  along  the  centerline  of  the  cylinder  normal  to  the 
laboratory  X\  axis  [i.e.,  X-axis  or  (230)]  is  shown. 
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Tables  6  and  7  list  maximum  local  cumulative  slip  (i.e.,  the  maximum  value  of  yk  at  any  location 
in  the  RDX  substrate)  at  an  indentation  depth  D  =  200  nm  when  the  respective  elastic  constants 
of  references  (/)  and  (2)  are  implemented.  Total  slip  y  listed  in  the  bottom  row  of  each  table  is 
not  necessarily  the  sum  of  all  yk  listed  in  a  given  column  because  the  location  in  the  substrate 
where  total  slip  is  maximum  does  not  necessarily  correspond  to  the  location  where  each  yk  is 
maximum.  Consistent  with  figure  13,  orientation  (001)  exhibits  slip  primarily  on  system  6. 
Significant  activity  of  multiple  slip  systems  is  evident  for  indentation  onto  (021)  and  (210)  planes. 
It  appears  that,  for  compatibility  reasons,  slip  system  6  is  active  for  indentation  on  (210)  because 
it  is  the  only  system  with  a  slip  direction  different  from  [100]  (see  figure  1  and  table  3).  Trends 
are  qualitatively  similar  regardless  of  choice  of  elastic  constants.  Cumulative  slip  magnitudes  are 
generally  slightly  larger  in  table  7  than  in  table  6  because  larger  stresses  are  attained  at  the  same 
depth  of  indentation  when  stiffer  elastic  constants  are  prescribed. 


Table  6.  Maximum  local  slip  yk  at  indentation  depth  of  200  nm  for  indentation  on 
(001),  (021),  and  (210)  planes,  elastic  constants  of  Haussuhl  (i). 


System  k 

Plane  (001) 

Plane  (021) 

Plane  (210) 

1 

0.01 

0.19 

0.40 

2 

0.01 

0.00 

0.40 

3 

0.08 

0.16 

0.29 

4 

0.08 

0.00 

0.29 

5 

0.00 

0.10 

0.51 

6 

0.31 

0.41 

0.16 

all  (y) 

0.31 

0.46 

0.52 

Table  7.  Maximum  local  slip  yk  at  indentation  depth  of  200  nm  for  indentation  on 
(001),  (021),  and  (210)  planes,  elastic  constants  of  Haycraft  et  al.  (2). 


System  k 

Plane  (001) 

Plane  (021) 

Plane  (210) 

1 

0.01 

0.19 

0.41 

2 

0.01 

0.00 

0.41 

3 

0.06 

0.13 

0.28 

4 

0.06 

0.00 

0.28 

5 

0.01 

0.10 

0.52 

6 

0.35 

0.51 

0.26 

all  (y) 

0.35 

0.56 

0.58 
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Loading-unloading  simulations  of  indentation  on  (001),  (021),  and  (210)  planes  of  RDX  are 
performed  to  enable  comparison  with  experimental  observations.  Force  versus  displacement 
predictions  are  shown  in  figure  28  for  a  single  load-unload  cycle  corresponding  to  loading  to 
D  «  200  nm  and  then  unloading  to  D  =  0.  During  much  of  the  loading  phase,  orientation  (001)  is 
slightly  less  stiff  than  orientations  (021)  and  (210).  However,  at  D  pa  200  nm,  P  is  largest  for 
orientation  (001).  Recall  that  orientation  (001)  demonstrates  a  lower  elastic  stiffness,  but  also 
less  slip  activity.  Hence  at  larger  indentation  depths,  increased  slip  activity  for  orientations  (021) 
and  (210)  lowers  their  effective  tangent  stiffness,  and  hence  P,  below  that  of  orientation  (001). 


D  [nm] 

Figure  28.  Indentation  force  vs.  applied  displacement  for  a  single  load-unload  cy¬ 
cle  on  each  of  (001),  (021),  and  (210)  planes. 


Hysteresis  is  also  substantially  greater  for  indentation  on  (021)  and  (210)  than  for  (001),  again 
demonstrating  less  slip  activity  in  the  latter.  Orientation  (210)  demonstrates  the  most  hysteresis 
(and  greatest  slip  activity)  and  the  largest  elastic  stiffness  for  much  of  the  unloading  phase. 

Residual  cumulative  slip  contours  for  indentation  to  a  depth  of  Dm  200  nm  onto  planes  (001), 
(021),  and  (210)  are  shown  in  figures  29,  30,  and  31.  In  each  case,  a  slice  along  the  centerline  of 
the  cylinder  normal  to  the  laboratory  X\  axis  (i.e.,  X-axis)  is  shown.  Contours  are  similar  to  those 
shown  in  figures  22,  23,  and  24,  but  with  greater  magnitudes  of  cumulative  slip  apparent  after 
unloading  on  (001)  and  (210).  Recall  from  equation  18  that  cumulative  slip  can  never  decrease  in 
time,  and  slip  in  forward  and  reverse  directions  on  a  given  system  both  contribute  positively; 
hence,  residual  slip  contours  can  demonstrate  a  larger  magnitude  as  a  result  of  any  forward  or 
reverse  slip  occurring  during  unloading.  Cumulative  slip  shown  in  figure  30  is  exceeded  by  that 
in  figure  23  since  in  the  former  the  depth  of  indentation  does  not  attain  200  nm  (figure  28). 
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Figure  29.  Residual  cumulative  total  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (001).  A  slice  along  the 
centerline  of  the  cylinder  normal  to  the  laboratory  X\  axis  [i.e.,  X-axis 
or  (100)]  is  shown. 


Figure  30.  Residual  cumulative  total  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (021).  A  slice  along  the 
centerline  of  the  cylinder  normal  to  the  laboratory  X\  axis  [i.e.,  X-axis 
or  (100)]  is  shown. 
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Figure  3 1 .  Residual  cumulative  total  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (210).  A  slice  along  the 
centerline  of  the  cylinder  normal  to  the  laboratory  X\  axis  [i.e.,  X-axis 
or  (230)]  is  shown. 
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Figures  32-42  show  total  cumulative  slip,  y,  and  cumulative  slip,  yk,  for  active  individual  slip 
systems  on  unloaded  surfaces  of  (001),  (021),  and  (210)  planes.  In  each  case,  the  maximum 
indentation  depth  prior  to  unloading  is  200  nm.  Non-circular  contours  are  consistent  with 
activity  of  fewer  than  five  geometrically  independent  systems  necessary  to  accommodate  an 
arbitrary  plastic  strain  field  (3). 


Figure  32.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D 
200  nm,  followed  by  unloading,  onto  plane  (001):  total  slip. 
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Figure  33.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (001):  slip  on  system  6 
(010)  [001]. 


Figure  34.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D 
200  nm,  followed  by  unloading,  onto  plane  (021):  total  slip. 
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Figure  35.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (021):  slip  on  system  1 
(021)  [100], 


Figure  36.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (021):  slip  on  system  3 
(011)[100], 
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Figure  37.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (021):  slip  on  system  6 
(010)  [001]. 


Figure  38.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D 
200  nm,  followed  by  unloading,  onto  plane  (210):  total  slip. 
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Figure  39.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (210):  slip  on  system  1 
(021)  [100], 


Figure  40.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (210):  slip  on  system  2 
(021)  [100], 
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Figure  41.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (210):  slip  on  system  5 
(010)[100], 


Figure  42.  Surface  residual  cumulative  slip  contours  for  indentation  to  depth  D  ~ 
200  nm,  followed  by  unloading,  onto  plane  (210):  slip  on  system  6 
(010)  [001]. 
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As  shown  in  figures  32  and  33,  surface  slip  activity  results  primarily  from  system  6  for 
indentation  on  (001).  Slight  contributions  to  the  circular  total  slip  trace  in  figure  32  are  due  to 
systems  3  and  4,  i.e.,  (01 1}  [100]  (not  shown).  As  shown  in  figures  34-37,  surface  slip  activity  is 
predominantly  from  system  6,  with  minor  contributions  from  systems  1  and  3,  for  indentation  on 
(021).  As  shown  in  figures  38-42,  surface  slip  activity  is  predominantly  from  system  5,  with 
contributions  from  systems  1,  2,  and  6,  for  indentation  on  (210).  Faint  contributions  from 
systems  3  and  4  (not  shown)  also  occur  for  indentation  on  (210). 

Active  slip  planes  during  the  loading  history  at  the  specimen  surface  as  predicted  by  the  present 
simulations  are  compared  with  those  deduced  from  experimental  surface  impressions  (3)  in  table 
8.  It  is  noted  that  in  simulations  of  indentation  on  (210),  systems  3  and  4,  i.e.,  (01 1}  [100], 
demonstrate  substantial  activity  within  the  bulk  material  (i.e.,  beneath  the  surface),  but  little 
cumulative  slip  at  the  surface.  Experiments  suggest  that  (01 1}  [100]  slip  traces  exist  for  indented 
(210)  surfaces  of  RDX  (3).  Though  not  presented  in  this  report,  predicted  residual  bulk  and 
surface  slip  contours  were  also  obtained  using  the  relatively  stiff  elastic  constants  of  Haycraft  et 
al.  (2).  Trends  in  results  for  slip  activity  were  usually  qualitatively  similar  regardless  of  choice  of 
elastic  constants  from  references  ( 1 )  or  (2). 

Table  8.  Significantly  active  slip  planes  at  the  specimen  surface  during  indentation 
of  (001),  (021),  and  (210)  planes. 


Indentation  surface 

Present  results 

Experiment  (2) 

(001) 

(010),{011} 

(010), {Oil}, {021} 

(021) 

(010), (Oil), (021) 

not  reported 

(210) 

(010), {Oil}, {021} 

(010),{011},(02l) 

4.  Conclusions 


A  nonlinear  anisotropic  elastic -plastic  model  has  been  developed  for  single  crystals  of  RDX.  The 
model  accounts  for  orthorhombic  elastic  constants,  pressure-dependent  compressibility,  and 
dislocation  glide  on  up  to  six  distinct  slip  systems  from  four  families  of  crystallographically 
equivalent  slip  systems. 

Numerical  simulations  of  spherical  indentation  on  (001),  (021),  and  (210)  planes  of  single 
crystals  have  been  performed.  Results  show  significant  influences  of  elastic  anisotropy  and 
elastic  nonlinearity  on  force-displacement  data.  Model  predictions  for  initial  elastic  response 
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using  constants  measured  with  resonant  ultrasound  spectroscopy  agree  with  experimental 
force-displacement  data  for  indentation  on  (001),  (021),  and  (210)  planes.  Model  predictions  for 
initial  elastic  response  using  constants  measured  with  Brillouin  scattering  are  in  reasonable 
agreement  with  experiments  for  indentation  on  (210),  but  are  stiffer  than  experiments  for 
indentation  on  (001)  and  (021).  Orientation  (001)  is  elastically  most  compliant  with  respect  to 
indentation  force,  in  agreement  with  experiments. 

Model  predictions  of  force  for  larger  indentation  depths,  wherein  local  plastic  slip  is  substantial, 
tend  to  exceed  experimental  values  regardless  of  which  set  of  elastic  constants  is  used.  The 
constant  strength  (i.e.,  perfectly  plastic)  slip  model  implemented  in  the  present  work  is  unable  to 
replicate  nearly  horizontal  steps  in  indentation  force  observed  in  experiments.  Such  steps  may 
correspond  to  discrete  slip  events  of  width  too  fine  to  be  captured  by  a  conventional  continuum 
slip  model.  Dependencies  of  shear  strength  and  shear  stiffness  on  slip  history  and  pressure  have 
been  omitted;  incorporation  of  such  physics,  for  example,  as  suggested  by  atomic  simulations, 
might  provide  improved  agreement.  Uncertainties  in  true  indenter  tip  geometry,  and  fractures 
observed  in  experiments  (at  the  surface)  or  not  observed  (subsurface)  could  also  contribute  to 
discrepancies  between  simulations  and  experiments. 

Critical  shear  strengths  associated  with  plastic  yielding  have  been  estimated  on  the  order  of  Gj 20, 
where  G  is  a  representative  elastic  shear  modulus.  Simulations  suggest  that  slip  planes  (010)  and 
{011}  contain  active  systems  for  indentation  on  (001),  with  slip  on  system  (010)[001]  dominating 
the  inelastic  response;  experimental  surface  observations  confirm  that  these  and  {021}  slip  planes 
may  also  be  active.  Simulations  suggest  that  slip  planes  (010)  and  {021},  and  to  a  lesser  extent 
{011},  are  active  at  the  specimen  surface  for  indentation  on  (210);  these  same  planes  have  also 
been  confirmed  as  active  in  experiments.  Simulations  suggest  that  planes  (010),  (Oil),  and  (021) 
contain  active  systems  for  indentation  on  (021);  particular  slip  planes  active  for  this  orientation 
have  not  been  reported  from  experimental  observations. 

The  present  simulations  suggest  that  system  (010)  [001]  provides  the  largest  contribution  to  the 
inelastic  material  response  (i.e.,  post-yield  force  versus  displacement  curve)  for  indentation  on 
(001)  and  (021)  planes,  while  five  systems  {021}[100],  {01 1 } [100] ,  and  (010) [100]  all  contribute 
to  inelastic  response  for  indentation  on  (210)  planes.  Plastic  deformation  and  hysteresis  are  more 
extensive  for  indentation  on  (021)  and  (210)  than  for  indentation  on  (001). 

Since  much  plastic  deformation  occurs  in  the  bulk  of  the  material,  and  since  different  slip 
mechanisms  may  be  prominent  at  the  surface  and  in  the  bulk,  these  results  offer  new  insight  into 
inelastic  mechanical  behavior  of  RDX  not  available  from  experimental  observations  of  residual 
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surface  topography  alone. 


In  summary,  the  model  developed  here,  when  used  with  elastic  constants  obtained  from  resonant 
ultrasonic  methods,  is  thought  to  provide  an  accurate  representation  of  the  nonlinear  anisotropic 
response  of  RDX  single  crystals  in  the  elastic  regime  of  indentation.  This  model  is  also  thought 
to  provide  a  qualitatively  realistic  depiction  of  activity  of  different  slip  systems  when  a  uniform 
and  constant  shear  strength  on  the  order  of  G/20  is  prescribed.  However,  refinements  of  the 
model  are  needed  to  address  any  reduction  in  stiffness  associated  with  discrete  or  highly  localized 
slip  events  and  cleavage  fractures  observed  at  larger  indentation  depths. 

Some  information  contained  in  this  report  appears,  in  abbreviated  form  with  far  fewer  figures,  in  a 
recently  submitted  publication  (37). 
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Appendix.  Rotated  Elastic  Constants 


Let  C^ljk]  be  the  fourth-order  matrix  of  elastic  constants  of  the  crystal  oriented  for  indentation 
into  crystallographic  plane  (z  jk).  Let  RWk)  be  the  corresponding  rotation  matrix.  Then 


(P ( ijk)  _  p(*A)  p( *' A)  p (ijk)n( ijk) r 
-'A BCD  ~  KAE  KBF  KCG  KDH  -'EFGH  ■ 


(A-l) 


where  C efgh  corresponds  to  the  isothermal  elastic  constants  of  table  2.  Let  C^j  denote  the 
rotated  elastic  constants  in  Voigt’s  notation. 


For  indentation  into  (001)  planes  of  RDX, 


Rm  =  [Rmh 


AB 


1  0  0 
0  1  0 
0  0  1 


(A-2) 


The  second-order  elastic  constants  of  Haussuhl  (/)  are 


_  [C(00t) 


24.56 


7.61 

18.85 


symm 
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17.33 
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The  second-order  elastic  constants  of  Haycraft  et  al.  (2)  are 
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24.49  8.16 

20.78 

symm 


0 

0 

0 

11.99 


0 

0 

0 

0 

2.72 


0 

0 

0 

0 

0 

7.68 


GPa. 
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For  indentation  into  (021)  planes  of  RDX, 


Rm  _  [^(021)1 


MS 


1  0  0 
0  0.475  -0.880 
0  0.880  0.475 


Rotated  second-order  elastic  constants  of  Haussuhl  (/)  are 


24.56  5.82 

16.78 


symm 


Rotated  second-order  elastic  constants  of  Haycraft  et  al.  (2)  are 


£(ijk)  —  p(021) 


36.48  1.18 

24.94 


symm 


7.09 

0.97 

0 

0 

6.13 

-0.27 

0 

0 

17.62 

0.90 

0 

0 

6.04 

0 

0 

6.26 

1.19 

4.70 

et  al. 

(2)  are 

0.98 

-0.15 

0 

0  ’ 

4.84 

2.96 

0 

0 

26.98 

-1.40 

0 

0 

8.67 

0 

0 

6.56 

2.07 

GPa. 


3.84 


GPa. 


(A- 5) 


(A-6) 


(A-7) 


For  indentation  into  (210)  planes  of  RDX, 


=  [*$0)] 


Rotated  second-order  elastic  constants  of  Haussuhl  (7)  are 

'20.15 


0.495 

-0.869 

0 

0 

0.869 

0.495 

cm  =  [cgj°) 


5.26 

7.72 

0 

1.16 

0 

17.33 

5.29 

0 

0.03 

0 

23.06 

0 

1.29 

0 

4.33 

0 

-0.47 

symm 

7.01 

0 

4.88 


GPa. 


(A- 8) 


(A- 9) 
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Rotated  second-order  elastic  constants  of  Haycraft  et  al.  (2)  are 


cm  =  [C(2W) 


22.17 


6.47 

20.78 


symm 


6.16 

0 

-0.54 

0 

2.95 

0 

-2.97 

0 

28.29 

0 

5.70 

0 

4.99 

0 

-3.99 

12.94 

0 

9.72 

A3 -direction,  C 

,(210)  r(021) 

33  ^  °33  ^ 

GPa. 


(A- 10) 


4001) 

'33 


for  either  set  of 


experimentally  measured  elastic  constants  for  RDX.  It  is  also  noted  that  C33  from  Brillouin 
scattering  reported  by  Haycraft  et  al.  (2)  always  exceeds  the  corresponding  value  from  resonant 
ultrasonic  measurements  reported  by  Haussuhl  (7). 


Rotation  matrices  and  elastic  coefficients  listed  in  equation  A-2-A-10  are  non-unique  since 
additional  rotation  about  the  normal  to  plane  (i jk )  would  not  change  the  crystallographic  plane 
normal  to  the  direction  of  indentation.  While  some  elastic  constants  would  change  under  such 
additional  rotation,  the  elastic  constant  C^Jk'  would  be  unaffected.  Furthermore,  for  spherical 
indentation  into  a  cylindrical  substrate,  the  total  indentation  force  P  would  be  unaffected  by  rigid 
rotation  of  the  crystal  lattice  about  A3. 

It  is  also  illustrative  to  consider  the  elastic  compliance  matrix  §  obeying  the  reciprocity  relations 

^ABEF^-EFCD  =  -  ( &AC&BD  +  5 AD§Bc)-  (A-l  1) 

Components  of  the  6  x  6  compliance  tensor  Sap  are  those  of  the  inverse  of  Cap .  Young’s  moduli 
are  defined  as 


£11  =  l/Sn  =  1/Smi,  £22  =  I/S22  =  l/§2222,  £33  =  I/S33  =  1/^3333-  (A- 12) 

Values  for  Young’s  modulus  for  orientations  corresponding  to  indentation  on  (001),  (021),  and 
(210)  planes,  with  A3  the  direction  of  indentation,  are  listed  in  table  A-l  for  each  set  of  elastic 
constants  given  in  table  2. 

Values  of  Eyik>  do  not  follow  trends  mentioned  for  those  of  C33  ,  and  trends  for  Young’s 
modulus  also  differ  between  different  sets  of  experimentally  measured  constants.  In  particular, 
for  elastic  constants  of  Hausshul  (7),  >  £33°^  >  £33”'  *  •  For  elastic  constants  of  Haycraft 
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et  al.  (2),  >  £33 101  >  £33° '  ‘  ■  It  is  remarked  that  E^a'  from  Haycraft  et  al.  (2)  always 

exceeds  the  corresponding  value  from  Haussuhl  (/)  for  a  =  1, 2, 3. 

Table  A-l.  Young’s  modulus  of  RDX  single  crystals  oriented  for  indentation  on 
planes  (001),  (021),  and  (210). 


Plane  [ref.(i)] 

Plane  [ref.(2)] 

Modulus 

(001) 

(021) 

(210) 

(001) 

(021) 

(210) 

En  [GPa] 

20.85 

20.85 

16.72 

36.40 

36.40 

19.12 

E22  [GPa] 

15.69 

14.05 

15.40 

21.28 

22.84 

18.04 

£33  [GPa] 

15.40 

14.19 

19.31 

18.04 

25.55 

23.61 
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